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ABSTRACT 

The atomic hydrogen gas (H i) disk in the outer region (beyond ~ 10 kpc from 
the centre) of Milky Way can provide valuable information about the structure 
of the dark matter halo. The recent 3-D thickness map of the outer H I disk from 
the all sky 21-cm line LAB survey, gives us a unique opportunity to investigate 
the structure of the dark matter halo of Milky Way in great detail. A striking 
feature of this new survey is the North-South asymmetry in the thickness map of 
the atomic hydrogen gas. Assuming vertical hydrostatic equilibrium under the 
total potential of the Galaxy, we derive the model thickness map of the H I gas. 
We show that simple axisymmetric halo models, such as softened isothermal halo 
(producing a flat rotation curve with V c ~ 220 kms -1 ) or any halo with density 
falling faster than the isothermal one, are not able to explain the observed ra- 
dial variation of the gas thickness. We also show that such axisymmetric halos 
along with different H I velocity dispersion in the two halves, cannot explain the 
observed asymmetry in the thickness map. Amongst the non-axisymmetric mod- 
els, it is shown that a purely lopsided (m = 1, first harmonic) dark matter halo 
with reasonable H I velocity dispersion fails to explain the North-South asym- 
metry satisfactorily. However, we show that by superposing a second harmonic 
(m = 2) out of phase onto a purely lopsided halo e.g. our best fit and more 
acceptable model A (with parameters e\ = 0.2, e\ = 0.18 and ghi = 8.5 kms -1 ) 
can provide an excellent fit to the observation and reproduce the North-South 
asymmetry naturally. The emerging picture of the asymmetric dark matter halo 
is supported by the ACDM halos formed in the cosmological N-body simulation. 



Subject headings: Galaxies: kinematics and dynamics - Galaxies: spiral - Galax- 
ies: structure - galaxies: ISM - galaxies: halos 
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1. Introduction 

Asymmetries are common in disk galaxies and are seen in lopsidedness (Baldwin et al. 
1980; Rix & Zaritsky 1995; Bournaud et al. 2005; Saha et al. 2007), in warps (Garcia- Ruiz 
et al. 2002; Sanchez-Saavedra et al. 2003; Saha & Jog 2006), in the rotation curves of 
receding side and approaching side (Manabe & Miyamoto 1975; Swaters et al. 1999), and 
in the distribution of the neutral hydrogen gas in galaxies (Richter & Sancisi 1994). These 
asymmetries in the dynamical phenomena can provide valuable information about the nature 
of the underlying dark matter potential. In fact, they have often been considered a reflection 
of the asymmetry in the dark matter distribution (Weinberg 1994; Jog 1997). However, 
one needs to be careful and systematically discard other possibilities such as dynamical 
instabilities as the origin of these asymmetries (e.g. Saha & Jog 2006). Our motivation in 
this paper is to check whether a non-axisymmetric dark matter halo can explain the recently 
measured asymmetry in the atomic hydrogen gas distribution in our Galaxy. 

In order to study the shape of the dark matter halo, it is best to look for tracers which 
lie mostly outside the main baryonic disk component. Examples of such tracers are the 
motion of satellites, neutral hydrogen (H i) gas in the outer region of the galaxy etc.. In 
particular, over the last two decades it has become well known that the H I gas layer flares 
significantly beyond the optical disk of our Galaxy (Kulkarni, Heiles & Blitz 1982; Knapp 
1987; Wouterloot et al. 1990; Diplas & Savage 1991; Merrifield 1992; Nakanishi & Sofue 2003) 
and the recent LAB (Leiden/Argentine/Bonn) survey (Kalberla et al. 2005) of Galactic H 
I reveals the most comprehensive, uniformly sampled flaring map of the H I gas extended 
out to a very large radius from the Galactic centre. Many external edge-on galaxies seem 
to show flaring in the thickness of the neutral hydrogen gas (Brinks & Burton 1984; Oiling 
1996; Matthews & Wood 2003). The most likely reason for the flaring is that the total 
gravitational force acting perpendicular to the disk plane decreases with radius while the 
velocity dispersion of H I is observed to be nearly constant (Lewis 1984). The contribution 
to the total perpendicular gravitational force comes mainly from the stellar disk, gas and the 
dark matter halo. Since the midplane density of the stars falls off rapidly compared to that 
of dark matter halo at large galactocentric radii (typically beyond the optical disk), the dark 
matter halo is expected to take over the major role in determining the vertical distribution 
of H I gas in the outer region. This makes the H I layer in the Galactic outskirts extremely 
sensitive to the distribution of dark matter and presently available neutral hydrogen gas 
(from the LAB survey) extended out to a very large radius from the Galactic centre provides 
us an unique opportunity to examine the detailed nature of the dark matter distribution in 
the Milky Way. 

This approach has been used to investigate the nature of dark matter halos by studying 
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the thickness of the neutral hydrogen gas in our Galaxy (Oiling & Merrifield 1998, 2000, 
2001, Narayan et al. 2005; Kalberla et al. 2007) and for M31 (Banerjee & Jog 2008) . Of 
these, recent work by Narayan et al. 2005 based on the Wouterloot et al. (1990) H I data 
shows that the observed H I flaring of Milky Way is best fit by a spherical dark matter halo 
with density falling faster than the isothermal halo. On the other hand, the most recent and 
very detailed work by Kalberla et al. (2007) , based on the LAB survey data, shows that 
the dark matter distribution in the Galaxy is rather complicated. They needed a massive 
extended dark matter halo, a self-gravitating dark matter disk, and a dark matter ring to 
explain the flaring in the H I gas in our Galaxy. Now, one of the most striking feature in the 
LAB survey is that the neutral hydrogen gas thickness map shows systematic North-South 
(hereafter N-S) asymmetry, where Galactic North refers to < ip < 180 ° and Galactic 
South by 180° < tp < 360°, in the Galctocentric cylindrical polar coordinate system (R, tp, 
z). This asymmetry was previously seen in previous HI maps of the Milky Way (Henderson, 
Jackson and Kerr 1982), but rarely remarked on. So it may be worth investigating if the 
Milky Way H I flaring can be explained by a non-axisymmetric dark matter halo. 

Such non-axisymmetries in the collisionless dark matter halo are probably not uncom- 
mon. In the current cosmological paradigm, the A cold dark matter (A CDM) halos are 
formed by dissipationless gravitational collapse of the material associated with the peaks of 
the primordial density fluctuation field and then grow via mergers and accretion in a highly 
nonlinear fashion. The resulting dynamical structures of these halos can be highly asym- 
metric. In fact, the A CDM halos formed in the recent Millennium Simulation (Springel et 
al. 2005) show asymmetry in their mass distribution (Gao & White 2006). 

In the present study, we derive, numerically, the flaring in the thickness of the neutral 
hydrogen gas using self-consistent model for the Galaxy (Narayan & Jog 2002; Narayan et 
al. 2005) including a non-axisymmetric dark matter halo and disk. We take into account the 
self-gravity of the gas. Our analysis produces a non-axisymmetric flaring curve for the gas 
and in this respect our study is probably different from all the previous studies which tried 
to derive H I flaring in spiral galaxies. Based on our analysis, we show that an elliptically 
perturbed (readers are referred to §2.2 for details) lopsided dark matter halo with density 
falling faster than that of an isothermal halo can explain the observed N-S asymmetry in 
the H I thickness map in the Galaxy. 

The paper is organized in the following order. §2 describes the formulation of the prob- 
lem and models of the dark matter halo. In §3 we present the definition of the thickness of 
the H I gas and asymmetry measurements. Method and input parameters are discussed in 
§4, while §5 describes the results and various possible models of dark matter halo. Com- 
parison of different models are done in §6. §7 describes a comparison with previous works. 
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Discussion and conclusions are in §8 and 9 respectively. 



2. Vertical dynamics in a Lopsided dark matter halo 



We studied the dynamics of an H I gas disk under vertical hydrostatic equilibrium in 
a generalized non-axisymmetric dark matter halo potential. The H I disk is gravitationally 
coupled to the stellar counterpart. Basically, we are going to study the equilibrium vertical 
structure of a two-component star-gas system in an asymmetric potential. We used the 
cylindrical polar coordinate system (R,tp,z) suitable for the disk geometry. Presuming that 
the H I gas is in hydrostatic equilibrium, the vertical dynamics of each component under the 
force field due to the surrounding dark matter halo can be described by coupling the Poisson 
equation and the equation of hydrostatic equilibrium along the normal to the mid-plane for 
each component. 

The Poisson equation for the system can be written as: 



where $j is the total potential due to the stars, H I gas and dark matter halo, p« with i=l 
to 2 denotes the mass density for the stellar and H I components, p h is the density of dark 
matter halo. 

Next, we made a comparative estimate of the Poisson equation's radial and azimuthal terms 
(which we denote as T R and T v respectively). We assumed that the potential of the Galaxy 
is asymmetric by a small quantity e p and the dominant asymmetry is in the form of a 
lopsidedness (i.e. the iso-potential contours are distributed according to cosip, corresponding 
torn = 1 azimuthal wavenumber). We then wrote the total potential, $ t , in the simplistic 
form (Rix & Zaritsky 1995; Jog 1997): 



where <&o(R, z) is the axisymmetric part of the total potential and <p p is the constant phase 
factor. Then the azimuthal term can be written as 




(1) 



<& t (R,ip,z) = %(R,z)(l + e p cos(y? - (p p )), 



(2) 




R 2 



e p cos(y? - (p p ), 



(3) 



and the radial term as: 
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1 dV 2 

Tr = R ~dR (1 + Cp C ° S(V? " ^ (4) 

where the rotation velocity, V c = a/ -R<9$o(-R, 2) / <9i? is determined at the mid-plane (z — 0). 
To proceed further with the comparison, we need the form of the combined axisymmetric 
potential §q(R,z). We assumed that the rotation curve in the outer region of the stellar 
disk is mainly dominated by the dark matter halo. 

Let V c ~ C R~ a , where the index a determines the shape of the rotation curve and 
C is some constant of proportionality. Note that a < i.e. a negative a would produce 
rising rotation curves usually found in dwarf galaxies. In normal spiral galaxies a > 0. For 
example, a = produces a flat rotation curve that corresponds to a screened isothermal 
dark matter halo (p — 1 in our notation, see eq.[12] in §2.2). In the other extreme, when 
a = 0.5 the rotation curve falls in a Keplerian fashion. The dark matter halos that produce 
an asymptotically Keplerian rotation curve are with indices 1.5 < p < 2 (For the present 
work, we restrict ourselves to p — 2). Therefore, we can say that the range of the indices 
in the rotation curves, < a < 1, roughly maps to the range 1 < p < 2 in the density 
distribution of the dark matter halos. Now clearly, for an asymptotically flat rotation curve 
the radial term is exactly equal to zero at the disk mid-plane (z — 0). Whereas, the azimuthal 
term is given by T v ~ log(i?)/i? 2 e p cos((/9 — (p p ). For a very slowly falling rotation curve, the 
azimuthal term dominates over the radial term and with little algebraic manipulations it can 
be shown: 



T R + [1 + J cos(v? - p p )]T R . (5) 

Then, the Poisson equation for a disk embedded in a non-axisymmetric dark matter halo 
with a slowly falling rotation curve («^0 and a C 1) can be written as: 

^ ~ 4nG Pi + P^j ~ [1 + J cos(^ - tp p )\T R . (6) 

Since T R is a negative quantity, the combined effect of the radial and azimuthal terms in the 
Poisson equation (Eq.[6]) is either to increase or decrease the vertical oscillation frequency of 
the disk depending upon the orientation of the dark matter halo and the values of a and e p . 
In other words, this term may either try to confine the gas more towards the disk mid-plane 
or help flaring. However, the contribution from the combined radial and azimuthal terms 
to the thickness of H I in the outer region is not significant compared to that due to the 
first term in eq.[6]. The second term on the r.h.s. of eq.[6] contributes to ~ 10% in the 
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H I thickness in the outer region for a very slowly falling rotaion curve. So in the zeroth 
order approximation, we can safely say that the thickness of the disk components is largely 
determined by the vertical pull due to the first term on the r.h.s. of the Poisson eq.[6]. 

The vertical hydrostatic equilibrium equation under the total potential for each component 
can be written as (Rohlfs 1977; Binney & Tremaine 1987): 

2 d\n Pi d$ t 

In the above equation a zi denotes the vertical velocity dispersion of the i th disk component 
in the problem. On combining the above equations (6, 7), the vertical equilibrium of each 
component in the disk under the dark matter halo potential is: 

' '■'(> r I \ ^cos[,- r> )]:/ ( ,,. (8) 



dz \ Zl dz 



-) = -47TG Pi + P^j + [1 + J cos(v? - <p p )\T R . 



2.1. An analytic approximation to the thickness 

An analytic approximation for the thickness can be obtained near the disk mid-plane by 
direct integration of the vertical equilibrium equation. This is a useful guide to the numerical 
integration of eq.[8]. 

On integrating eq.[7] of the vertical equilibrium equation and using the boundary condition 
at z=0 

Pi(R,ip,z) = p l mid {R^,Q) 
, where p l mid (R, f, 0) is the mid-plane volume density of each component, we get: 



a,i In 



Zl 



+ $t(R,<p,z)-$ t (R,<P,0) = 0. (9) 



Near the disk mid-plane the potential along the vertical direction will not be very different 
than that of the disk mid-plane. Using Taylor's expansion along the vertical direction and 
the vertical equilibrium of each component, it can be shown that: 

Pi (R,<p,z) = PIM^O) x e ~^#^>. (10) 
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Hence, the density distribution under vertical equilibrium follows a gaussian near the disk 
mid-plane with non-axisymmetric thickness. Here Hi(R,(p) = a zi /iy(R,<f) and the vertical 
frequency v is given by 



u(R,ip) 



^G[J2pln ld (R^) + Ph(R,v,o) 



;n) 



In deriving the above form, we used eq. [6] with T R = 0, since the thickness of the gas is 
largely determined by the vertical pull of the matter towards the disk mid-plane. The reader 
should bear in mind that the above analytical formula for the thickness (Hj) is not valid for 
high z. 



2.2. Models of Non-axisymmetric Dark Matter Halo 

We consider here for the first time a seven-parameter dark matter halo model for de- 
scribing the neutral hydrogen gas thickness in the Milky Way. The density profile of the 
non-axisymmetric dark matter halo is given by: 

Ph(R,<P,z) = r P0(q l hp , (12) 



A 2 

1 + 



R 2 M 

where p is the central mass density of the halo, R c is the core radius, q determines the 
oblateness of the halo and p is the index determining the nature of the density profile. The 
index p — 1 denotes a softened isothermal dark matter halo with density falling as Rr 2 at 
large radii. The name 'softened isothermal halo' is derived from the singular isothermal halo 
by adding a suitable core radius; these softened halos are also called screened or pseudo- 
isothermal halos. Since the mass M(R) within a radius R is proportional to R, the p — 1 
halo produces an asymptotically flat rotation curve. For p — 1.5, the density ph oc R~ 3 at 
large radii resembling the NFW halo profile (Navarro et al. 1996). The mass M(R) within 
a radius R of such a p — 1.5 halo is proportional to logi? and goes to infinity as R goes to 
infinity but much more gradually than the p — 1 halo. Whereas p = 2 denotes a perfect 
ellipsoid dark matter halo with density falling like R~ 4 at large radii leading to essentially a 
finite mass halo. So the asymptotic rotation curve would be like a Keplerian one. The mass 
profile of the p = 2 axisymmetric dark matter halo can be written as: 
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M 2 (R) ~ 27T Po R i c q tan'^R/Rc) 



R/R c 



(13) 



1 + R 2 /R 2 C 



In the above eq. [12] 



Al = {Rl + z 2 /q 2 )- R v = R(l - e\ cos (<p 



(fi h ) + 4 cos2(ip -cp h )), 



(14) 



where represents the surface of the concentric ellipsoid with lopsided (m = 1) distribution 
superposed with second harmonics (m = 2). The degree of lopsidedness in the dark matter 
distribution is e\ and the degree of second harmonic is determined by t\. iph denotes the 
phase of the asymmetric halo with respect to the Galactic axes. It is very hard to say without 
a detailed stability analysis which harmonics (m = 1 or m = 2) will dominate in the dark 
matter halo. 

Since the higher harmonics (m > 2 or 3) perturbations will be associated with smaller 
spatial scales compared to the first or second harmonics, in a collisionless dark matter halo 
such small scale higher harmonics are most likely to be Landau damped and large scale 
perturbations would be weakly damped (Weinberg 1994). Our expectation is that the dark 
matter halos are likely to be dominated by the first two harmonics: a lopsided (m = 1) per- 
turbation and an elliptical (m = 2) perturbation. The effects of a lopsided halo perturbation 
and m = 2 component in the halo perturbation onto the disk dynamics was studied by Jog 
(1999, 2000). Note that with e^=0, we end up with a purely lopsided dark matter halo. 
With all €h = 0, the dark matter halo becomes the usual four-parameter axisymmetric halo 
used in previous studies (Narayan et. al. 2005; de Zeeuw & Pfenniger 1988; Becquaert & 
Combes 1997). The parameter iph contains important information as to how the dark matter 
halo is oriented with respect to the Galactic axes. It is important to note that the first two 
harmonics in the dark matter halo are out of phase. Writing the density distribution of the 
halo in the above form results in freedom to investigate a wide variety of dark matter halo 
potentials, while simpler to use than the triaxial ones. 



By solving the coupled Poisson equation and the vertical hydrostatic equilibrium equa- 
tion (eq.[8]), we obtain the volume density of the atomic hydrogen gas p 2 (R,(p,z). To solve 
eq. [8] we used the mid-plane volume density and its dispersion of each component as inputs 
(see §4). We considered an exponential surface density distribution for the stellar disk with 
a mild lopsidedness into it. 



3. Thickness of the neutral hydrogen gas 
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V(R, if) = S e-^(1 + e s cos(^ - <p h )\ 



(15) 



where S and are the central surface density and scale length of the stellar mass dis- 
tribution, respectively. The value of e s , in principle, can be determined self-consistently by 
calculating the response of the stellar disk to the non-axisymmetric dark matter halo consid- 
ered in §2.2 but the task is beyond the scope of this paper and not necessary for obtaining a 
correct zeroth order model. Instead, we use reasonable numerical values of e s and calculate 
its effect on the H I thickness in the region of our interest (R > 16 kpc). We have checked 
that the gas thickness remains almost unchanged as the value of e s is reduced from 0.15 to 
0.0 (corresponding to an axisymmetric stellar disk). We have also checked if e s alone can 
account for the observed systematic N-S asymmetry (prominent beyond ~ 16 kpc) in the H 
I thickness map. We found that even unrealistically high values of e s alone can not account 
for such high asymmetry in the H I flaring. Note that the contribution of the axisymmetric 
stellar disk itself almost drops to zero beyond ~ 16 kpc compared to that due to the dark 
matter halo. So in the present calculation, we considered only the axisymmetric stellar disk 
in deriving the H I thickness beyond ~ 16 kpc. One important assumption in solving eq. [8] is 
that the velocity dispersion remains constant (isothermal approximation) along the vertical 
direction. 

The thickness d(R, </?) of the gas is determined by using the second moment of the 
volume density distribution and is given by the following relation: 



d(R, ip) gives the radial variation of the thickness along a particular azimuthal direction (</?) 
in the disk. Using this thickness map, we can examine the degree of asymmetry in the gas 
thickness distribution. 

We note that the thickness of the gas layer in this case is different from the method 
used by Levine et al. (2006a) but gives qualitatively similar result. Kalberla et al. (2007) 
calculate the gas scale height yet another way, and get result which differ both from the 
second moment of the distribution (eq. [16]) and from Levine et al. (2006a). All this says is 
that because of the exclusion of certain areas (especially 90 < <p < 110°), and the possible 
effects of optical depth, different methods of calculating the scale height can give rise to 
different numerical values. In fact, it was demonstrated by Bahcall (1984) that the scale- 
height of a gaussian stellar distribution is roughly twice that of an expoential distribution. 
The important point, is that all methods show values of scale-height that are ~ 2 — 2.5 times 



d 2 (R^) 



!°° O0 z 2 p 2 (R^ ) z)dz 
f™ OQ p 2 (R,(p,z)dz 



(16) 
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larger in the northern part of the Galaxy (1 = — 180°) than the south (/ = 180 — 360°) 
and it is this difference that we trying to model and explain. 

Let A(<f k ) be the area under the thickness curve along a particular azimuthal angle ((p k ): 

A(ip k )= / d(R,cp k )dR (17) 

In the above equation R m i n and R max represent the initial and the final radius respectively 
in the available observation of the gas thickness. The degree of asymmetry in the thickness 
distribution, presuming that we are considering a smooth curve, can then be described by: 

vM - %)T%^)' (18) 

The range of r\ is < r] < 1, with 77=0 denoting a symmetric distribution. On the other 
hand i]=l denotes a highly asymmetric gas/star distribution leading to a one-sided flaring 
in the galaxy. It is unlikely that rj ~ 1, because this would mean that the velocity dispersion 
of the gas in one half of the galaxy is almost zero (cold component) compared to the other 
half. This may lead to an instability in the disc. Considering (p k = 90°, the above relation 
would produce the numerical value of the observed N-S asymmetry in the thickness map. 



4. Method and Input parameters 

The present study focuses on deriving the thickness map of the neutral hydrogen gas (H 
I) in the very outer region (R > 16 kpc) of the Galactic disk. Since of the various baryonic 
components the stars and H I dominate the mass in this region, we neglect the effect of 
molecular hydrogen gas (H 2 ) on the thickness distribution of H I. Eq.(8), which describes 
the zeroth order vertical equilibrium under a non-axisymmetric dark matter distribution, 
represents two coupled differential equations for the two disk components: H I and stars. 
The vertical density distribution for each component, responding to the total potential due 
to the disk and the dark matter halo, was solved numerically as an initial value problem 
using the fourth order Runge-Kutta method of integration (Press et al. 1994). The details 
of this method are presented in Narayan & Jog (2002) and Narayan et al. (2005). Because 
of the underlying non-axisymmetry, eq.(8) is solved for each azimuthal direction along the 
Galactocentric radius. Along a particular azimuth, at a particular radius, we calculate the 
second moment of the vertical density distribution for each component according to eq.(16) 
and call it the thickness of that component. Repeating this procedure for regular intervals 
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in the azimuth along the galactocentric radius give us the thickness map of the neutral 
hydrogen gas layer in the Galaxy. 

The primary input parameters needed for the method are the mid-plane volume density 
and the vertical velocity dispersion for each disk component. The H I mid-plane volume 
density is obtained from the LAB survey (Kalberla et al. 2005; Levine et al. 2006a). For 
the stellar disk, we used its surface density according to eq.(15) for which we need to know 
the central surface density and scale length. Note that eq. [15] represents a lopsided stellar 
disk and as discussed in §3, the value of e s does not effect much the H I thickness in the 
very outer region of the Galaxy. So we considered effectively an axisymmetric stellar disk 
in the region of our interest. Using the following measured/inferred quantities: the stellar 
surface density at the solar region E , the disk scale length R d and the distance of sun from 
the Galactic center R , we can derive the central surface density S and infact the surface 
density at any radius. We use £ = 45 M pc~ 2 which is consistent with 48 ± 9 M pc~ 2 
obtained by Kuijken & Gilmore (1991) and 52 ± 13 M pc~ 2 obtained by Flynn & Fuchs 
(1994) for the total surface density, after the gas density is subtracted. We use the IAU 
recommended value for R & (=8.5 kpc) and this was also used to determine the observed 
thickness map for H I gas (Levine et al. 2006a). The scale length was set equal to 3.2 kpc 
(Mera et al. 1998) in accordance with the recent determinations of smaller disk scale-length 
for our Galaxy. 



4.1. Stellar and H I velocity dispersion 

The stellar vertical dispersion was derived from observation of radial dispersion by Lewis 
& Freeman (1989) and then using the assumption that the ratio of the vertical to radial 
velocity dispersion is equal to 0.45 at all radii in the Galaxy, equal to its observed value in 
the solar neighbourhood as obtained from the analysis of the Hipparcos data (Dehnen & 
Binney 1998, Mignard 2000). 

The H I velocity dispersion {(Thi) has been observed to be nearly constant with radius 
at approximately 9±1 km s _1 (Spitzer 1978; Malhotra 1995) in the inner Galaxy (out to the 
solar circle). Beyond the solar circle, however, the dispersion has not yet been measured. A 
study of 200 external galaxies (Lewis 1984) shows that the observed dispersion has a very 
narrow range, about 8±1 km s -1 , consistent with observations of our Galaxy. Sicking (1997) 
showed that in two external galaxies, dispersion decreases slowly out to the outer edge of the 
H I layer. In a number of other galaxies, the velocity dispersion decreases and then stabilizes 
at a constant value of 7±1 km s _1 (Shostak & van der Kruit 1984; Dickey 1996; Kamphuis 
1993). This decrease in velocity dispersion is perhaps due to the lesser number density of 



- 12 - 



supernovae (SN) in the outer region (McKee & Ostriker 1977). On the other hand, recent 
work by Dib et al. (2007) shows that the SN do not affect the observed gas velocity dispersion 
in the galactic outskirt. A large fraction of the observed velocity dispersion is non-thermal (or 
turbulent) in origin and the supernovae could be the major source for this turbulent nature 
of the H I velocity dispersion. Note that the thermal contribution accounts for only about 1 
km s _1 (Spitzer 1978). In any case, the H I velocity dispersion is as crucial as the dark matter 
distribution is in determining the H I thickness map in the outer Galaxy. Unfortunately, in 
the absence of any direct measurement of the H I velocity dispersion (a hi) beyond the solar 
circle, we use am as a model parameter in the fitting problem. We construct various model 
based on the H I velocity dispersion values and models of dark matter distribution (Eq.[12]). 



4.2. Model rotation curves 

We model the rotation curve for the Galaxy using a bulge, an exponential disk for the 
stars and gas and a dark matter halo. Modelling the observed rotation curve is a non-trivial 
task given that there are various kind of uncertainties and turns out to be a non-linear 
regression problem in a multidimensional parameter space. Naturally, such modelling would 
suffer from uniqueness problem. So we aim here to reproduce the main features in the 
observed rotation curve and try to generate such rotation curves for which the circular speed 
lies in the range determined by the relation O = (27±2.5)-R kms -1 due to Kerr & Lynden- 
Bell (1986). Since we use the IAU recommended value for R Q = 8.5 kpc, the above range 
implies V c = 230 ±21 kms^ 1 at the solar radius. Keeping these constraints in mind, we 
proceed to derive the model rotation curves in the following way. We assume that the disk 
and the bulge are aligned with the symmetry axis of the dark matter halo. We also assume 
the virial equilibrium in the Galaxy. Then the square of the total circular velocity in the 
disk mid-plane (z = 0) can be written as: 

v c = Vl lge + V s ] ars + V g 2 as + V* mh , (19) 

where we adopt a Plummer-Kuzmin bulge model to derive the bulge contribution to the 
rotation curve. The density profile of the bulge is given by the following formula (Binney & 
Tremaine 1987): 

where R b is the bulge scale-length and M b is the total bulge mass 



(20) 

. We have used same 
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values of these two parameters in all of our model rotation curves e.g. Rb = 2.5 kpc and 
Mb = 2.8 x 1O 1O M (see Blum 1995). The inclusion of this simple spherical bulge reproduces 
a reasonable looking rotation curve for the Galaxy. However, the contribution of the bulge 
to the H I thickness in the region of interest (R > 16 kpc) is almost negligible, the bulge 
contributes to ~ 5% to the overall gas thickness. 

The disk contribution to the circular speed is derived using the potential for the expo- 
nential mass distribution (S(i?) = £ e where values of the parameters are mentioned 
in §4) and similarly for the various dark matter halo models for which the parameters are 
mentioned in the appropriate places below. The gas contribution to the total circular speed 
is derived from the following exponential distribution for the H I beyond 14 kpc (Levine et 
al. 2006a): 

E HI (R) = 4.5 X e [-(*-14fcpc)/4.3fcpc] 

Below 14 kpc, we use a constant surface density for the H I gas (as observed) to derive 
its contribution to the rotation curve. Although there is a notable difference of our model 
rotation curves with the observed one (Brand & Blitz 1993), the slowly falling rotation curves 
due to the p = 2 halo model follow the trend found in the Milky Way's rotation curve in the 
very recent analysis of Xue et al. (2008) based on the SDSS data. 

5. Hi flaring and nature of dark matter halo 

There is a clear North-South asymmetry in the thickness map of the H I gas in our 
Galaxy. It is not obvious what would have caused such asymmetry in the H I gas distribution. 
What is the underlying nature of this asymmetry? It remains to be shown whether this is 
purely a gas dynamical effect or reflects some gravitational effect. Below we describe a 
step by step analysis of the cause of this asymmetry, which gradually reveals the nature of 
dark matter halo in our Galaxy. We used averaged flaring data over the north (0 < < 
180°) and south (180° < 4> < 360°) respectively excluding about 15° region about the 
Sun-Galactic centre line. The vertical volume density distribution of the gas is derived using 
the rotation curve due to Brand & Blitz (1993). The thickness of the gas is then derived 
by taking the second moment of the density distribution after the thickness filter has been 
applied. In this respect, the thickness measurement of the gas is different from the Levine 
et al. (2006a) who uses the tail integration method and it is also different from Kalberla et 
al. (2007) who uses half width at half maximum (HWHM) for the half-thickness of the gas 
later. In all our analyses below, we use the same thickness map for the gas, although it is 
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understood that a different rotation curve would indeed produce different thickness map for 
the gas. However, on doing an error analysis we find that the error in calculating the distance 
due to an error in the rotation velocity is small (~ 10% at large galactocentric distances) to 
produce an appreciable change in the observed thickness of the gas. The readers are referred 
to §8 for a discussion on the dependence of the thickness data on the rotation curve. The 
observed N-S asymmetry denotes the average asymmetry in the thickness map of the H I 
gas and according to eq.[18], its given by i] obs = 0.262. 

5.1. Models of axisymmetric dark matter halo 

We first considered a simple axisymmetric dark matter halo model (with all = in 
eq.[14]) for the Galaxy and used different velocity dispersions for the H I in the two halves 
to see if the observed N-S asymmetry could be reproduced. 

5.1.1. p — 1, Softened Isothermal Halo 

The softened isothermal dark matter halo produces naturally the asymptotically flat 
rotation curve in a spiral galaxy. It has been shown previously based on the Wouterloot 
et al. (1990) data that a p = 1 softened isothermal halo of any shape (oblate or prolate) 
cannot explain the observation (Narayan et al. 2005). In the present study, we found the 
same trend, so we considered a nearly spherical halo to begin with for further investigation. 
We adopt the parameters of the p — 1 isothermal dark matter halo from the mass model 
of our Galaxy based on microlensing observation by Mera et al. (1998). We found that 
the p = 1 halo with core radius R c — 5 kpc, central density po=0.035 M Q pc -3 , which gives 
rise to an asymptotically flat rotation curve with terminal velocity of 220 km s" 1 (Brand 
& Blitz 1993) (see bottom right of Fig. 1), and with <Jhi = $ km s _1 cannot explain the 
present observation (LAB data). The fact that p — 1 softened isothermal halo can't explain 
the flaring in the gas thickness distribution has already been verified by Narayan et al. 
(2005) in the case of Wouterloot et al. (1990) data and more recently by Kalberla et al. 
(2007) based on the LAB survey data. Now, there are two aspects of the H I thickness 
map derived from the LAB survey data: the radial variation of the gas thickness in the 
North and the South and the prominent asymmetry between the North and the South as 
mentioned already. An axisymmetric p — 1 halo with the above mentioned parameters 
(giving rise to V c = 220 km s _1 ) and <jhi = 9 km s _1 appears to be quite massive and hence 
difficult to reproduce the present observation. By inspection of Fig. 1 (top left), it is easy 
to check that there has to be an order of magnitude decrease in the dark matter mid-plane 
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density to reproduce the observed gas thickness at 30 kpc. Moreover, the striking difference 
between the two slopes makes it clear that even a linearly increasing velocity dispersion with 
radial distance (although unphysical) cannot fit the observation. Obviously, an axisymmetric 
p = 1 halo alone can not explain the other important aspect of the thickness map, namely, 
the observed N-S asymmetry. This fact leads us to explore another possibility of a model of 
an axisymmetric p — 1 halo accompanied by a non-axisymmetric distribution of H I velocity 
dispersion. It can be explained in a simple way that a different H I velocity dispersion on 
both sides of the hemispheres is not going to improve the situation either. In a naive theory, 
the thickness of H I gas can be written as hni oc cfhi / y/Gp m id- Using this simple formula, it 
is easy to see that in order to increase the thickness by a factor of 2 at a particular radius, 
one needs to raise the dispersion by a factor of ~ 2 (because p m id remains the same). At 
R = 30 kpc, the thickness of H I gas in the Northern hemisphere is roughly a factor of 2 more 
than that in the Southern hemisphere and to explain this asymmetry based on purely gas 
dynamical effect and the axisymmetric p — 1 isothermal halo (giving rise to a flat rotation 
curve), one needs to have a hi ~ 20 km s _1 which is unlikely according to the standard 
models of ISM of our Galaxy. 



5.1.2. p=1.5, an NFW type halo 

At large distances (R 3> R c ) from the centre, p = 1.5 halo resembles an NFW profile 
(Navarro et al. 1996) for the dark matter halo. In our model of vertical equilibrium, this 
seems to be preferable compared to the p = 1 softened isothermal dark matter halo. Because 
of the density falling faster than the p = 1 isothermal halo, the resulting rotation curve also 
starts falling beyond about 10 kpc. We use eq. [2.96b] of Binney & Tremaine (1987) to 
generate the rotation curve of the p = 1.5 halo. The axisymmetric p=1.5 halo with R c = 8 
kpc and p =0-025 M pc -3 produces a reasonable rotation curve (see bottom right of Fig. 1). 
At Rq the rotation velocity is 223 km s _1 and at 2R & , V c =216 km s -1 . At 25 kpc, the 
difference in rotation velocities between the p = 1 and p = 1.5 halo is approximately 20 km 
s -1 . In the top right panel of Fig. 1, we show the half-thickness of H I gas due to the p = 1.5 
halo model considered here. The axisymmetric halo model with <jhi=9.2 km s~ x fits well 
the observation in the southern part beyond about 15 kpc. However, the same model does 
not fit the thickness curve at all in the northern halves. The solid line shows the curve with 
H I velocity dispersion ~ 12.2 km s _1 and yet does not give a good fit. An increase in the H 
I velocity dispersion beyond 12 km s _1 would force the model curve to intersect the observed 
one only at a single point in the North. This again demonstrates that the N-S asymmetry in 
the thickness map of the H I is probably not due to a gas dynamical effect, rather it arises 
due to some gravitational effect. 
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Fig. 1. — Half thickness of neutral hydrogen gas (first three panels) in the Milky Way ana- 
rotation curve models (bottom right) for different dark matter halos. Top left: Solid line 
is the model fitted to the North and dashed line to the South for the p — 1 axisymmetric 
halo giving rise to a flat rotation velocity of 220 km s^ 1 at large radii. Top right: Model 
half thickness due to axisymmetric p — 1.5 dark matter halo with different a hi in the two 
halves. Bottom left: Model half thickness due to axisymmetric p = 2 dark matter halo with 
different am in the two halves 
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5.1.3. p = 2, a Perfect Ellipsoidal Halo 

Previous studies by Narayan et al. (2005) have investigated the axisymmetric p = 2 
halo in considerable detail to explain the H I thickness in the Galaxy. Their p = 2 halo 
provides a good fit to the Wouterloot et al. (1990) flaring data. In the present study, we also 
considered an axisymmetric p=2 halo and see if different H I velocity dispersions can explain 
the N-S asymmetry. We consider a core radius R c — 9.4 kpc and central density p =0.035 
M pc~ 3 which produces a reasonable rotation curve (see Fig. 1) within the uncertainties in 
the observation out to 30 kpc where the difference between the p — 1 (flat rotation curve) 
and p = 2 rotation curves is ~ 25 km s _1 . Note that the rotation curves for both p — 1.5 
and p = 2 halos are slowly falling with radial distance. Recent analysis of the kinematics 
of a large number of Blue Horizontal-Branch halo stars from the SDSS database by Xue et 
al. (2008) show that the rotation curve of our Milky Way is actually falling slowly with the 
Galactocentric radius. So this direct observational analysis supports the falling trend in the 
rotation curve of Milky Way. Infact, the choice of the halo parameters for the p=2 case 
makes the rotation curve more realistic as observed. At 25 kpc, the circular speed due to our 
p = 2 halo differs from the flat rotation curve only by ~ 15 km s _1 which is well within the 
observed error bars. Beyond about 15 kpc, with linearly decreasing H I velocity dispersion 
from 9.2 (at 10 kpc) to 8.0 (at 30 kpc) km s -1 , the p = 2 halo considered here gives a good 
fit to the observed data in the southern halves. However, even with ani—Ylfl or 14.2 km 
s -1 , the p = 2 halo does not fit well the observation in the north. With 14.2 km s" 1 , the 
model curve over-estimates the observation between 13 - 23 kpc and underestimates beyond 
24 kpc. On the other hand with 12.2 km s -1 , the model under-estimates the observation 
beyond 20 kpc. The solid curve in Fig. 1 (Bottom left panel) is with a constant a HI =13.2 
km s -1 in the Northern halves. Compared to p — 1 and p = 1.5, an axisymmetric p = 2 halo 
with different H I velocity dispersion comes closer to the observation. We have also tried to 
use a linearly decreasing H I velocity dispersion from 13.2 or 12.2 km s -1 , but no good fit 
found. It becomes quite clear that the observed H I flaring in the North cannot be explained 
with any physically meaningful variation of the H I velocity dispersion. 



5.2. Non-axisymmetric dark matter distribution 

Based on the above three axisymmetric cases we have investigated so far, we conclude 
that an asymmetric state of ISM combined with axisymmetric dark matter potential is not 
able to explain the observed asymmetry in the thickness map indicating that the observed 
asymmetry is probably not due to a gas dynamical effect. The large scale asymmetry in the 
north-south thickness map is instead likely to be gravitational in origin. Beyond 16 kpc, 
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the stellar contribution is not significant compared to the dark matter, which implies that 
the N-S asymmetry is not likely due to the asymmetric stellar disk and also the observed 
asymmetry in the stellar disk is not very high in the Galaxy. This gives us a room to explore 
whether the N-S asymmetry originates due to an asymmetric dark matter halo. As explained 
in §1, such asymmetric dark matter halos are not uncommon even in the cosmo logical N- 
body simulations. By inspecting the observed H I thickness map, we guess that the N-S 
asymmetry is primarily lopsided (m = 1) in nature. A lopsided dark matter halo has been 
used to provide an explanation for the observed lopsidedness in the underlying stellar disk 
and for the asymmetry in the rotation curves of disk galaxies (Jog 1997, 2002). Below we 
test our hypothesis that the Milky Way's dark matter halo is lopsided. We show next that 
the observed nature of the H I distribution demands that the lopsided dark matter halo to 
be oriented with phase angle <ph = 270°. In the present case, fa is the angle between the 
direction in which the dark matter distribution is elongated more on one side compared to 
the other and the Galactocentric coordinate axes. So the dark matter density maximum 
of the lopsided halo is along the southern direction (0 = 270°) and the density minimum 
is along the northern direction (0 = 90°). We test p = 1.5 and p = 2 halos as the likely 
candidates in our further investigation; p — 1 is not used because it fails to produce the 
observed flaring even in the south. 
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Fig. 2. — Half thickness of neutral hydrogen gas in Milky Way. Solid line is the model fitted 
to the North and dashed line to the South. The value of Hi velocity dispersion a m is 9.2 
kms -1 . 



-19- 



5.2.1. p=1.5 lopsided halo 

Having oriented the dark matter halo with phase <ph = 270°, we found that an ax- 
isymmetric (e^=0) p = 1.5 halo with H I velocity dispersion 9.2 km s _1 fits quite well the 
observation in the southern halves (see Fig. 1). However, a purely lopsided configuration 
(e\ = in eq.[14]) with reasonable parameters can not explain the observed asymmetry in 
the thickness map of H I gas. The parameter ranges we have tried are a hi = 7—10 km 
s _1 and e\ = 0.02 — 0.4. We found that no combination of these parameters could provide 
a good fit to the observed data. Further, we tried to fit the combined model of the dark 
matter halo (given in eq. [12]) with lopsidedness as the dominant component. Again, for 
the above mentioned parameter range no good fit to the observed data was found. Fig. 2 
shows a configuration of the dark matter halo in which both the first (m = 1) and second 
(m = 2) harmonics are equal in strength and out of phase with each other. Of course, such 
configuration of p — 1.5 halo with reasonable H I velocity dispersion (9.2 km s _1 ) provides 
better fit to the observation than the axisymmetric configuration of p — 1.5 halo with high 

5.2.2. p=2 lopsided halo 

Since the axisymmetric p = 2 halo already seemed promising compared to others, we 
carried out a detailed analysis of the p = 2 lopsided halo here. We used the basic parameters 
like core radius and central density of the axisymmetric p = 2 configuration (see §5.1.3). In 
order to find a best fit model of the Galaxy under a purely lopsided halo, we make a 2D 
grid of two independent and free parameters (<Jhi, €-\) spanning a large dynamical range i.e. 
[cr hi ^ 4H( 7 " 12 ) km s_1 , (°- 05 - °- 40 )] for tnis halo model. We found that with a H i < 9.0 
km s" 1 , a purely lopsided p = 2 halo is not sufficient to explain the observation. With some 
more exploration, we found that a purely lopsided p=2 halo with e\ = 0.17 and a hi = H 
km s _1 does fit the observation quite well (see Fig. 3). We call this configuration model p2L. 
The rotation curve for this model is shown in the bottom panel of Fig. 3. The difference 
between the rotation velocities in the North and the South is ~ 20 km s" 1 for this model. 
The density contours of the dark matter halo are shown in the right panel of Fig. 3. Note 
that the H I velocity dispersion for this model p2L is fairly high. So we further explored 
the parameter space consisting of e\, e\ and a hi for the p = 2 dark matter halo to find the 
best possible model to explain the observed radial variation of the thickness map and its N-S 
asymmetry. 

Based on the different values of H I velocity dispersion and different values of the e\ and 
e\, we construct three models (model A, model B and model C) which give a very good fit 
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Fig. 3. — Left: Half thickness of neutral hydrogen gas in Milky Way, model p2L. Solid 
line is the model fitted in the North and dashed line to the South. See Table 1 for the 
model parameters. Right: Density contours of the p = 2 lopsided dark matter halo, dashed 
lines are for the axisymmetric perfect ellipsoidal (p — 2) halo. Solid lines are for the purely 
lopsided p = 2 dark matter halo in case of model p2L. R c is the core radius of the dark 
matter halo. Contour levels are 0.05 xp for the inner most one and decreasing by a factor 
of 2 outwardly. Hi velocity dispersion a hi is 11.0 kms -1 for this model. 



to the averaged observation in both halves of the Galaxy and thereby explain the observed 
asymmetry. 

Model A 

We assumed the velocity dispersion of H I, a hi = 8.5 km s _1 , to be constant out to 
30 kpc for simplicity. The value of the velocity dispersion is more reasonable one and it is 
very closed to what has been used by Kalberla et al. (2007) in their best fit model. Once 
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a hi is fixed, we are essentially left with again a 2D grid of free parameters e\ and e\ for 
the dark matter halo. The parameter range explored is [e\, e^] = [(0.05 - 0.40), (0.05 - 0.40)]. 
We found that the inclusion of the second harmonic component in the dark matter halo 
dramatically improved the quality of the fit compared to all the previously explored cases. 
The best fit model parameters, for the above mentioned a hi, are 6^=0.20 and e^=0.18. 
Such configuration of the dark matter halo which is indeed lopsided, reproduced the N-S 
asymmetry quite well in the Galaxy (see Fig. 4). The rotation curves for the model in the 
two halves of the Galaxy are shown in the right panel of Fig. 4. The circular speeds at the 
solar radius are quite comparable to the p — 1 flat rotation curve within the uncertainties. 
At 25 kpc, the north and south rotation curves differ from each other by ~ 19 km -1 , the 
value in the north being ~ 209 km" 1 . The density contours of the model dark matter halo 
are shown in the bottom panel of Fig. 4. Because of the presence of the second harmonics 
(m = 2), the contours are elongated in the East- West direction of the Galaxy too. It is quite 
clear now that a simple model that consists only of a lopsided dark matter halo does not 
reproduce the observed asymmetry in the thickness map of the atomic hydrogen gas in our 
Milky Way. The lopsided dark matter halo with some amount of second harmonic (m — 2) 
superposed onto it appears to give the best fit to the data compared to the purely lopsided 
dark matter halo. We call the emerging picture of dark matter halo as elliptically perturbed 
lopsided dark matter halo. 

Model B 

Here we consider the velocity dispersion of H I, am = 9 km s _1 slightly higher compared 
to model A and again flat out to 30 kpc. The increased velocity dispersion improves the 
fit by a very small amount below 16 kpc. With the same halo parameters as in model A, 
the fit is not good beyond 16 kpc. So we reduced the strength of the second harmonic (e 2 h ) 
and when the second harmonic is decreased to ~ 2/3 of the first harmonic (ejj we again 
recovered a good fit to the data. The resulting model halo with parameters e\ = 0.2 and 
e\ = 0.14 reproduced the N-S asymmetry in the Galaxy quite well (see Fig. 5). The rotation 
curves for the model in the two halves of the Galaxy are shown in the right panel of Fig. 5. 
At 25 kpc, the rotation speed in the north is ~ 191 km -1 and in the south it is ~ 212 km -1 . 
It appears that the rotation curve in the south is more closer to the flat rotation curve (with 
V c = 220 km -1 ). The density contours of the model dark matter halo are shown in the 
bottom panel of Fig. 5. Due to the presence of a small component of second harmonic, the 
contours are elongated by a small amount compared to model A in the East- West direction 
of the Galaxy. This dark matter halo can be considered as a dominantly lopsided halo. This 
model again confirmed our earlier findings from model A. 
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Fig. 4. — Top left: Half thickness of neutral hydrogen gas in Milky Way, model A. Solid 
line is the model fitted in the North and dashed line to the South. Beyond ~ 16 kpc the 
fit is quite good. The model parameters are given in Table 1. Top right: Rotation curve of 
the Galaxy for model A. Bottom: Density contours of the dark matter halo of model A. 
dashed lines are for the axisymmetric perfect ellipsoidal (p=2) halo. Solid lines are for the 
elliptically perturbed lopsided p=2 dark matter halo with second harmonics (m = 2) being 
~ the first (m = 1) one (see table 1). R c is the core radius of the dark matter halo. Contour 
levels are 0.05 xp f° r the inner most one and decreasing by a factor of 2 outwardly, a hi 
used for this model is 8.5 kms -1 . 
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Fig. 5. — Top left: Half thickness of neutral hydrogen gas in Milky Way model B. The 
solid line is the model fitted in the North and dashed line in the South. Beyond ~ 16 kpc 
the fit is quite good. The model parameters are given in Table 1. Top right: Rotation curve 
of the Galaxy for model B. Bottom: Density contours of the dark matter halo of model 
B. dashed lines are for the axisymmetric perfect ellipsoidal (p = 2) halo. Solid lines are for 
the asymmetric lopsided p = 2 dark matter halo with the second harmonic (m = 2)~ 50% 
of the first (m = 1) one (see table 1). R c is the core radius of the dark matter halo. Contour 
levels are 0.05 xp for the inner most one and decreasing by a factor of 2 outwardly. a HI 
used for this model is 9.0 kms -1 . 
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Fig. 6. — Half thickness of neutral hydrogen gas in Milky Way, model C. Solid line is the 
model fitted in the North and dashed line to the South. The H I velocity dispersion is quite 
high in the radial range 10-16 kpc. See Table 1 for the dark matter halo parameters. 

We constructed this model primarily to see if we can explain the observation also in the 
region 10 < R < 16 kpc. We found that the mid-plane density distributions of H I and stars 
are almost same in the radial range 10 < R < 16 kpc (see Fig. 8) and interestingly enough 
the half-thickness of H I in this radial range is almost the same as that of stars. Another 
important fact to be noted is that the H I half-thickness is roughly constant in this radial 
range (10 < R < 16 kpc). Now, we know that the thickness of H I ~ cthi/ \f(Gpmid), where 
p mic i is the total mid-plane volume density. So we vary the H I velocity dispersion in the 
following manner in this radial range: 



From 16 kpc (~ 5 R d ) onwards a HI =8.5 km s" 1 and stays flat out to 30 kpc. We consider 
R a = 4370 pc because a similar value in the case of stellar disk can produce a flat thickness 
(Lewis & Freeman 1989) in the Galaxy. 

With the H I velocity dispersion varied in the above fashion (eq.[21]) and t\ = 0.2 and 
e|=0.18 (same as model A), we got an excellent fit to the observation (see Fig. 6). Model 
C provides the best fit to the observed data and reproduces the observed asymmetry in the 
thickness map of H I. In terms of halo asymmetry parameters, model C is exactly equal to 
the model A and so are the rotation curves. The density contours are the same as in Fig. 4. 
So model C confirms that our best fit dark matter halo model and again it is a lopsided dark 
matter halo with some elliptical perturbation. 



a HI = SS.ee^/^kms- 1 for 10 < R < 16 kpc . 
= S.Skms" 1 R > 16 kpc 



(21) 
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6. Comparison of different models 

6.1. Measurement of asymmetry 

Given the radial variation of the H I thickness in different azimuths ((f)), we can quantify 
the degree of observed asymmetry by rj (see eq.[18]) in the two halves of the Galaxy. We 
used this parameter r\ as a measure of the underlying asymmetry in the thickness distribution 
and to differentiate between the above models which seem to reproduce the observation quite 
well. The value of i]^ is computed by assuming = 90 ° as discussed at the end of §3. Just 
by comparing the observed asymmetry (?7 & s ) and the model predicted asymmetry (r) mo d e i), 
we find that model p2L is the best match to the observation (see Table 1) but the velocity 
dispersion used in this model is quite high. Our next best model is model C based on the 
comparison of asymmetry and this model also has a very high velocity dispersion between 10 
- 16 kpc. On the other hand model A is probably more reasonable because the H I velocity 
dispersion is not very high. The H I dispersion of model A is almost the same as that used 
recently by Kalberla et al. (2007); their solution is close to a single component model like 
we have considered with an effective constant velocity dispersion of 8.3 km s _1 . In all the 
four models of p = 2 dark matter halo one thing is common and it is the parameter e\ ~ 0.2, 
the degree of lopsidedness. So the different models are essentially built up around the same 
basic configuration, a lopsided halo. To make the picture more clear, these models, being 
a function of two apparently independent parameters (a hi and e\), can be thought of as 
a family of models of lopsided dark matter halos. From Table 1, amongst the first three 
models, it is clear that as e\ increases, a hi decreases to make the fit better. One interesting 
fact about these family of models is that a better fit to the observed data demands that the 
second harmonic be out of phase with the first one. 



Table 1: Asymmetry and the models of the p = 2 lopsided dark matter halos 



Model 


0~HI 
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Vobs 


Vmodel 




(km s _1 ) 


M pc- 3 


kpc 








M pc- 2 






p2L 


11.0 


0.035 


9.4 


0.95 


0.17 


0.0 


72.0 


0.262 


0.261 


A 


8.5 


0.035 


9.4 


0.95 


0.20 


0.18 


80.0 


0.262 


0.265 


B 


9.0 


0.035 


9.4 


0.95 


0.20 


0.14 


78.0 


0.262 


0.271 


C 


Eq.[21] 


0.035 


9.4 


0.95 


0.20 


0.18 


81.0 


0.262 


0.260 



In Table 1 i] obs is obtained from the observed average thickness of H I on both the halves. 
Vmodei is from the fitted curves. 
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6.2. Constraints from the surface mass density at Solar radius 

Apart from the constraint on the rotation curve due to Kerr & Lynden-Bell (1986), we 
provide here our estimations of total surface density at the solar radius for various models and 
compare them with other measurements from the literature. The value of the total surface 
density at the solar radius put a strong and important constraint on any mass model of our 
Galaxy. Based on K dwarfs, Kuijken & Gilmore (1991) provide us the total surface density 
within 1.1 kpc from the Galactic midplane E L i = 71 ±6 M pc~ 2 . By analysing HIPPARCOS 
K giants, Holmberg & Flynn (2004) measure Ei.i = 74 ± 6 M pc~ 2 . By modeling the H I 
thickness in the Galaxy, Kalberla et al. (2007) estimate Ei.i = 79 ± 2 M pc~ 2 . Our best 
fit models are in good agreement with all the previous measurements. For example, for 
the model p2L Ei.i = 72 M Q pc~ 2 ; for the model A, Ei.i = 80 M pc -2 ; for the model B, 
Ei.i = 78 M Q pc -2 and for model C, Ei.i = 81 M pc~ 2 . We used the total baryon mass 
surface density at the solar radius to be E star+9(1S = 50 M pc -2 which is within the limit 
(48 ± 8 M pc~ 2 at R ) provided by Kuijken & Gilmore (1989). Within 800 pc from the 
Galactic plane Holmberg & Flynn (2004) reports a surface density of E .s = 65 ± 6 M pc -2 
and from Kalberla et al. (2007) E .8 = 66.9 ± 2 M pc~ 2 . In this context, our model p2L 
gives E .8 = 65 M pc -2 and from model A, we get E .8 = 71 M pc~ 2 . In Fig. 7, we plot the 
total surface mass density (derived from modeling of the vertical hydrostatic equilibrium) of 
the Galaxy (including stars, gas and dark matter) within 1.1 kpc from the Galactic midplane 
as function of radius. It shows clearly that the surface density is higher in the Southern part 
of the Galaxy. In any case, our results on the surface density measurements in the solar 
neighbourhood are in close agreement with previously quoted values or within the quoted 
error bars. 
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Fig. 7. — The total surface mass density in the Galaxy within 1.1 kpc from the Galactic 
midplane. It includes the stars, H I and the dark matter distribution from our model p2L. At 
2R , the difference in the total surface densities between the North and South is ~ 4M pc~ 2 . 
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6.3. Total mass of the Galaxy 

Starting from the cosmological point view to studying galactic dynamics and especially 
in Galactic astronomy it is very important to have a good knowledge of the total mass of 
Milky Way. One of the complete analysis using most of the usual methods e.g. satellite 
kinematics, very high velocity stars, Local Group and Galactic rotation curve, Kochanek 
(1996) estimated the total mass of the Galaxy to be ~ 4.9 ± 1.1 x 10 11 M© within 50 kpc 
from the center of the Galaxy. Previous to this study, Lin, Jones & Klemola (1995) found 
the mass of the Milky Way to be 5.5 ± 1.0 x 10 11 M Q within 100 kpc from the centre. By 
modeling the rotation curve, Dehnen & Binney (1998) provided the Milky Way's mass within 
100kpctobe7±2.5xlO n M . From the uvby - (3 survey of high velocity stars, Garcia Cole 
et al. (1999) suggest the mass of Milky Way to be 3.2 — 3.6 x 10 11 M within 50 kpc from the 
Galactic centre. By using the kinematic information for Galactic satellites and halo objects, 
Sakamoto et al. (2003) estimated the mass of the Galaxy within 50 kpc (excluding Leo I) 
to be 5.0 — 5.5 x 10 11 M Q . In contrast to these studies, Kalberla et al. (2007) based on the 
H I thickness modeling in the Galaxy found the total mass of the Milky Way disk itself to 
be 2.9 ±0.1 x 10 11 M Q . The most recent analysis based on the Blue Horizontal-Branch halo 
stars from the SDSS, Xue et al. (2008) estimate the mass of our Galaxy to be 4.0 ± 0.7 x 10 11 
M within 60 kpc from the Galactic centre. Based on the modeling of H I thickness, we find 
the total mass of Milky Way to be ~ 3.3 x 10 11 M within 100 kpc from the center and 
within 50 kpc it is about 3.1 x 10 11 M . The mass of stellar disk is ~ 4.2 x 10 10 M n and the 
halo mass within 100 kpc is ~ 2.53 x 10 11 M . The difference between the halo masses in the 
North and South is ~ 2 — 3 x 10 10 M in our models of asymmetric halos. It is true that the 
total mass of the Galaxy saturates beyond about 100 kpc because the density profile of our 
p=2 dark matter falls faster (oc R~ 4 at R » R c ) than the the isothermal one. But within 
100 kpc, our estimates are in good agreement with most of the previous mass estimates of 
Milky Way. 

7. Comparison with Previous Work 

The most recent work which has extensively used the LAB survey H I data are Levine 
et al (2006a) and Kalberla et al (2007). Of these two, Kalberla et al. (2007) has studied 
in considerable detail the behaviour of the H I thickness in the Galaxy and found that in 
order to explain the observation they needed beside a massive dark matter halo, a dark 
matter disk and a dark matter ring. Such configurations of the dark matter for our Galaxy 
is troublesome for the CDM paradigm. 

Here, we would like to point out few similarities and differences of our method with the 
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previous ones. First of all, one of our best-fit models, namely model A, uses a constant H I 
velocity dispersion of 8.5 km s _1 which is close to what Kalberla et al. (2007) used for the 
their model namely 8.3 km s _1 . To make it more clear, although we do not use explicitly the 
different phases of ISM (e.g. CNM and WNM) to model the H I thickness map, the state of 
the ISM is roughly the same in our model as it is in Kalberla et al. (2007). In contrast to 
Kalberla et al. (2007), we consider non-axisymmetric models of dark matter halo and also 
a mild lopsided stellar disk. From Fig. 8, it is clear that the gas self-gravity is comparable 
to the stellar one and we take into account the full self-gravity of the H I gas in deriving the 
thickness map of the gas. 

There is a marked difference in the averaged flaring curves in the North and the South 
derived in our paper with that in Kalberla et al. (2007) which excludes a region 90 < (p < 
110 ° in the North showing very high flaring in the H I thickness map. On the other hand, if 
this region is included in the averaging process, the averaged observed gas thickness would 
increase by ~ 20% and if the H I velocity dispersion remains constant, one would expect a 
less massive dark matter halo from a simple minded calculation. Similarly, if we exclude this 
region the thickness in the North reduces to ~ 3.4 kpc. 

Another important fact about the LAB survey data is the pronounced N-S asymmetry in 
the gas thickness distribution. Certainly, an axisymmetric model of the Galaxy can not 
reproduce such asymmetry in the data. In this context, recent work by Sanchez- Salcedo et 
al. (2008) have also shown that MOND provides a reasonably good fit to the azimuthaly 
averaged flaring in the H I gas using the same LAB survey data. Now for an axisymmetric 
baryon distribution, the MOND potential would also be axisymmetric, because the differen- 
tial operator acting on the potential in the Poisson equation is rotationally invariant. Hence, 
it would probably be very hard to reproduce the observed N-S asymmetry in the H I thick- 
ness map. Whereas our family of lopsided dark matter halos naturally explain the observed 
North-South asymmetry. 

The density distribution of the dark matter halo (namely the p = 2) in our model is similar 
to that obtained by Narayan et al. (2005) (namely a p = 2 halo as the best fit model) based 
on the Wouterloot et al. (1990) H I flaring data. The axisymmetric version of our model is 
similar to that used by Narayan et al. (2005). 

8. Discussion 
(i) Self-gravity of the atomic hydrogen gas 

In the radial range from 10 - 16 kpc, we find that azimuthally averaged mid-plane volume 
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Fig. 8. — Mid-plane volume density of stars and gas in the two halves of the Galaxy. Density 
distributions are shown for the model A. For other models they are almost the same. Note 
that H I gas is comparable to the stars. N= North and S= South 



density of the atomic hydrogen gas is comparable to that of the stars (see Fig. 8) and beyond 
this radial range gas is in fact dominating over the stars. This suggests that the self-gravity 
of the H I gas is quite important in determining its thickness and thereby the nature of dark 
matter halo. By self-gravity of the gas we mean that the gas is held by its own gravity. In 
other words, we include the gas density in the Poisson equation to derive its contribution to 
the total potential of the Galaxy. Without any self-gravity, the gas would move like a test 
particle under the potential of the Galaxy. So it would be interesting to check the change 
brought about by excluding the H I self-gravity on its vertical scale- height. The difference is 
not negligible, it is seen to be about 10-15% within the optical disk (R < AR d ). This could be 
because these regions are dominated by the stellar disk and by the dark matter respectively. 
For the region 4 < R/Rd < 6, the difference is substantial (~20 - 30%) suggesting that in 
this range, the gas gravity is very important in negotiating the hydrostatic equilibrium for 
the H I layer. Thus neglecting it may lead to a serious overestimate of the H I scaleheight 
in general at all radii in the outer Galaxy and to explain the observed gas scale-height one 
may need to invoke a heavier dark matter halo. 

(ii) Molecular hydrogen gas 

In the present work, we are mostly concerned with the very outer region of the Galactic 
disk especially beyond 16 kpc. The particular reason for this is that the thickness curve 
of the atomic hydrogen gas in the northern hemisphere of the Galaxy deviates noticeably 
from the southern one beyond this region. The molecular hydrogen gas extends upto about 
17 kpc in the Galaxy (Wouterloot et al. 1990) and beyond around this region there is little 
data. Since our calculation of gas thickness is local (in the sense that we do not consider 
the gas gravity in global sense) we have neglected the effect of molecular hydrogen on the 
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H I thickness throughout our work. In fact, the absence of molecular hydrogen gas beyond 
about 16 kpc makes it more convenient in disentangling the effect of dark matter halo on 
the H I gas. 

(iii) Effect of Galactic constants 

We have built the stellar disk model based on the IAU-recommended values for the 
galactic constants - R — 8.5 kpc and O = 220 km s _1 . It would definitely be interesting 
and also worthwhile to know how the results for the halo density profile vary with the assumed 
galactic constants. For example, Oiling & Merrifield (2001) find the effect of varying these 
constants on the inferred axis ratio of the halo. Unfortunately, all the observational inputs 
for our model, like the H I and H 2 surface densities, H I scale-height and the stellar velocity 
dispersion, are based on the IAU-recommended galactic constants and rescaling them for 
different values of the constants is beyond the scope of this paper. 

(iv) Rotation curve and gas thickness data 

Different rotation curves produce different degrees of flaring because of the dependence 
of R on through equation 1 of Levine et al. (2006a). The rotation curves used in the present 
analysis differ by a small amount < 10% from the flat rotation curve, Q(R) = 220 km s _1 , 
used by Levine et al. (2006a). An error analysis shows that R changes approximately by 
~ 0/00, where is the circular speed at the solar position. Thus, dR/dQ, the sensitivity 
of an error in the distance to an error in the circular speed produces an error in the distance 
of a parcel of gas from the center of about 10% at 30 kpc. Because all values of Galactic 
latitude are small at large R, this translates to an error in the thickness of the gas layer by 
no more that 10%, and then only at the largest distances. 

Although the differences in the rotation curves used by us and by Levine et al. (2006a) 
implies that the two analyses are not exactly commensurate, we are trying to explain a factor 
of 2 increase in the thickness of the gas layer of the Milky Way from one hemisphere to the 
other, an order of magnitude larger than the maximum 10% effect caused by differences in 
the rotation curves. Note that the factor of 2 change in the scale height occurs at all radii 
beyond about R = 17 kpc and the effect on AR will be smaller at smaller radii because A0 
is also smaller. Therefore, while there is a 10% effect at the largest distances, the most it 
will do is to have a 10% effect on the scaling of the model. 

It is worth mentioning at this point that we do not solve the vertical hydrostatic equi- 
librium considering epicylic orbit correction for the gas. For this one needs to solve the Jeans 
equation and Poisson equation self-consistently, which is considerably more complex than 
what we do here, and in any event probably results in only a small correction. Our primary 
aim is to understand the nature of the observed asymmetry and build a first order model to 
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explain it. 

(v) Effect of Galactic warp 

The vertical structure of the Milky Way's disk is fairly complicated. The warp in the 
H I disk is very asymmetric; in the northern hemisphere the disk mid-plane rises to a height 
of about 4 kpc, while in the south it goes down to about 1 kpc and then again comes back 
to the undisturbed mid-plane (Kerr 1957; Burton 1988; Levine et al. 2006a). Since the 
warp is a global feature in the Galaxy, it will have its own self-gravity to affect the vertical 
oscillations in the disk. The effective vertical oscillation frequency in the disk can be written 
as v e ff = \jv 2 d + ^ warp , where v d is the vertical frequency of the unperturbed disk and v warp is 
the extra vertical frequency due to the enhanced self-gravity of the warp. Because the warp 
is global in nature, v marv is an integral quantity (see eq. [4b] in Saha & Jog 2006). However, 
an analytic form for v warp can be written in a local sense via the WKB analysis and it is 
~ a/27tGS^[aJ| (Binney & Tremaine 1987), where \k\ is the wavenumber for the warp. Since 
we are solving the vertical hydrostatic equilibrium locally, the contribution from the global 
warp (\k\ — > 0) becomes insignificant. In order to account for the warp in the gas thickness, 
one has to formulate the vertical hydrostatic equilibrium as an integral problem which we 
plan for the future. 

(vi) Galactic spiral structure 

The recent work by Levine et al. (2006b), based on the 21 cm LAB Galactic H I survey, 
reveals a multi-armed spiral structure of our Galaxy. Their study shows a good correlation 
between the positions of the spiral arms and the H I thickness. This is an expected behaviour 
because the presence of the spiral arms would cause enhanced gravity and scattering of the H 
I clouds would not change the vertical velocity dispersion appreciably, leading to a decrement 
in the H I thickness. However, we do not expect the thickness of H I to change appreciably 
because the strength of the perturbed surface density does not vary strongly as a function of 
the Galactocentric radius along an arm or even from arm to arm (Levine et al. 2006b). In 
the present paper, under the zeroth order approximation, we worked with the average data in 
the northern and southern halves of the Galaxy respectively. This is certainly an incomplete 
modeling of the data and in future we expect to construct a self-consistent formulation of 
the problem to include the spiral structure of the Galaxy. 

(vii) Uniqueness of the Lopsided halo model 

In astronomy, it is quite a hard job to prove directly the uniqueness of a proposed model 
which fits a given set of observational data well. One obvious way of approach is to compare 
different theoretical models against the given data set. Given the various uncertainties in 
the observation and our incomplete understanding of the physics of our Galaxy, such a 
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comparison may even lead to degeneracies and bayesian analysis would probably be the 
appropriate way out to lifting such degeneracies or discard other obvious models based on 
physical grounds. Before we go into discussing other possible theoretical models, we would 
like to remind the reader some basic facts just for a recap. We are dealing with the outer 
region of the Galaxy and especially beyond R > 16 kpc which is ~ 5 disk scale length and 
this is also the typical size of the stellar disk in a disk galaxy. So the influence of the stellar 
disk on the vertical distribution of the gas is almost negligible compared to that due to the 
dark matter halo. 

Certainly, the internal disk instabilities alone (being normally weak, except the bar 
instability) are unlikely to produce the observed asymmetry in the gas thickness distribution. 

One of the possible alternative models is an off-centered dark matter halo with respect 
to the Milky Way's disk. Such a model of an off-centered axisymmetric halo has been used 
previously by Levine & Sparke (1998) to generate lopsidedness in the disk. As pointed out 
by these authors that this method is effective when the disk lies within the core radius 
(almost constant density region of the halo) of the halo and as result is more efficient in 
dwarf galaxies rather than in luminous galaxies. In this context, one of our halo model 
namely p = 2 halo whose core radius is ~ 9 kpc as compared to the disk scale lenth ~ 3 kpc 
could have been a possible case for investigation. However, note that the N-S asymmetry in 
the gas distribution begins beyond 16 kpc which is roughly ~ 2 x R c of the halo and at this 
radius, the disk is no more in a nearly constant density region of the halo, making it harder 
to maintain the lopsidedness. In any case, this is an interesting possibility to be investigated 
in a future problem. 

Other possible model is an asymmetric gas accretion onto the disk as proposed by Bour- 
naud et al. (2005) to reproduce lopsidedness observed in the galactic disk. To produce strong 
lopsidedness as observed, one needs the accretion of the cold gas through the cosmological 
filaments to be highly asymmetric and in reality it is not clear if the gas accretion is asym- 
metric to such a degree. On the other hand, it makes sense to think in this direction because 
the H I distribution looks more disturbed in the North rather than in the South. At this 
point, it is worthwhile to mention that the gas actually contributes very little to the rota- 
tion curve (mostly dominated by the dark matter); so an axisymmetric p = 1 halo with an 
asymmetric gas accretion again may not be the good candidate for the present observation. 
However, an axisymmetric p = 2 halo with an asymmetric gas accretion could have been a 
possible candidate, but its beyond the scope of the present paper to examine such a model 
in considerable detail. 

On the other hand, our lopsided halo models are more natural to occur in cosmological 
scenario. Tidal interactions or large scale tidal harassment or major mergers of neighbouring 
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dark matter halos are more likely to produce large scale perturbations in the dark matter 
halo. And since the CDM halos are collisionless object the survival of such global perturba- 
tions in the halo is less of a problem. 

9. Conclusions 

We have analyzed both axisymmetric and non-axisymmetric configurations of dark mat- 
ter halos with a consistent picture of the ISM to explain the observed nature of the North- 
South asymmetry in the thickness distribution of the H I gas. Below we summarize the main 
results: 

With a model of the ISM that has reasonable values of the gas velocity dispersion, an 
isothermal dark matter halo producing a flat rotation curve with V c =220 km s" 1 cannot 
produce the observed flaring in the H I gas thickness. 

We show that the nature of the systematic North-South asymmetry in the H I thickness map 
is gravitational in nature. An axisymmetric dark matter halo with different values of the 
H I velocity dispersion in the two halves of the Galaxy can not reproduce this asymmetry. 
The observed asymmetry in the thickness map of neutral hydrogen gas is apparently not the 
result of purely gas dynamical effects. 

We show that a purely p — 1.5 or p — 2 lopsided dark matter halo also cannot explain 
the observed North-South asymmetry. For a purely p = 2 lopsided halo, the H I velocity 
dispersion has to be unreasonably large to come close to the observation and even then, the 
fit is not very good. 

Finally, we come up with a configuration of the dark matter halo in which some amount 
of second harmonic (m = 2) is superposed out of phase onto a purely p = 2 lopsided halo. 
For the best fit models (A & C), the values of lopsidedness and elliptical perturbation are 
e\ = 0.2 and e\ = 0.18 respectively. We call such a halo an elliptically perturbed lopsided 
dark matter halo which can explain the observed North-South asymmetry. Basically, the 
emerging picture of the dark matter halo of the Milky Way is dominantly lopsided in nature. 
In such a halo, the density falls off faster than the p — 1 isothermal dark matter halo. The 
emerging model of the asymmetric dark matter halo is supported by the halos formed in the 
recent cosmological N-body simulation. 
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